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Abstract 

We present the analysis of the phase space geometry of 2 — ► 3 reaction for the 
general case of nonzero and unequal particle masses. Its purpose is to elaborate an 
alternative approach to the problem of integration over phase space which does not 
exploit the Monte Carlo principle. The fast and effective algorithm of integration 
based on Gauss method is developed for treating 1-dimensional distributions in 
two-particle invariant variables. The algorithm is characterized by significantly 
improved accuracy and it can meet requirements of interactive processing. 
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1 Introduction 

Analysis of kinematics and integration over phase space of a reaction are basic problems 
which are common to both experimental and theoretical high energy physics. 

In the theoretical field the integration over phase space, for example, is needed to 
make theoretical predictions confronting the experimental measurements and to provide 
fitting of experimental data. Besides, a lot of pure theoretical investigations such as 
statistical models of nuclear reactions, calculation of unitarity corrections, etc., rely upon 
the integration over phase space of final state. 

At all stages of the modern experiments from designing the experimental device to 
determination of its characteristics and, finally, data treatment, various kinematic tools, 
including integration packages, are used for simulation of events, determination of accep- 
tances and making corrections to the measured distributions. 

The Monte-Carlo integration principle built into the most elaborate GEANT3 
system for modeling experimental devices provides subroutines for integration purposes 
which are universal in the sense that any number of final particles can be treated. The 
time-of-run and precision characteristics of the Monte-Carlo based programs are not so 
crucial here since, for example, determination of acceptances and other device properties 
is the direct problem. At the same time the majority of theoretical applications are in the 
field of inverse problems — to find out physical parameters some kernel built of theoretical 
model and device characteristics (angles, momentum cuts, etc.) must be inverted when 
applied to the experimental results. It is well-known that the precision requirements for 
determination of the kernel properties which obviously include the procedure of integration 
over a domain of phase space are much more stringent than in the case of the direct 
problem. 

In contrast with experimental applications where the integration is usually being per- 
formed with empty phase space the theoretical ones operate with complicated model 
amplitude which itself needs considerable time for processing. The time-of-run proper- 
ties of the integration routines become especially vital for creating an interactive utilities 
for data bases and routines for data treatment (like phase shift analysis package SAID 

0)- 

Therefore, development of an alternative to Monte-Carlo methods of phase space 
integration even at the price of loosing universality in respect to the number of particles 
seems to be important. 

Up to now alternative to Monte-Carlo methods are used only in the simplest case of 2 
— > 2 reactions (the phase space is effectively 1-dimensional and, therefore, any integration 
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method, in fact, applies well in this occasion). Appearing from time to time papers like 
||. [| , heavily rely on simplification of phase space geometry due to equal masses or specific 
amplitude of the reaction. 

Meanwhile the case of the 2 — > 3 reactions which importance many times was stressed 
since the work || by Chew and Low admits almost equally complete elaboration as in 
the 2 — > 2 case. The ground for such a statement might be found in the many-years-long 
investigations the results of which are summarized in the monographs || p]]. Of course, 
the full-scale realization of a classical integration scheme for all choices of variables of 
2 — > 3 reactions and at any kind of kinematical constraints will require an approach of 
artificial intelligence and computers of extraordinary power. 

The main goal of the present research is the demonstration of the principal feasibility 
of the classical integration approach. This will be done by developing the integration 
algorithm for treating total cross sections and 1-dimensional distributions in two-particle 
invariant variables for the general case of unequal particle masses of 2 — > 3 reaction. This 
problem requires extended analysis of the geometry of the phase space of the reaction. 

Fortunately, almost all the necessary examinations and hints might be found in various 
particular chapters of the books 110]'- one needs to gather carefully crucial conclusions 
and put considerations into purposive line of algorithm. Certainly, all the key relations 
must be expanded in exact analytic form to provide solid programming. The presentation 
of principle answers is an incidental goal of the paper. 

Another principle goal is to find a way beyond the standard integration algorithms. 
The extremely effective Gauss integration method is proved to be quite feasible in the 
considered case. 

The paper is organized as follows. Sect. 2 introduces basic notations and the starting 
form of the phase space integral. The central Sect. 3 is devoted to analysis of necessary 
elements of geometry of the phase space. To avoid burst of formulae most of explicit 
expressions are concentrated in Appendix. The iterated form of the phase space integral 
is introduced in Sect. 4 where the most important properties of the integrated expressions 
are investigated to motivate the applicability of Gauss method. The specific characteristics 
of the FORTRAN implementation of the results of the previous sections are considered in 
Sect. 5. The perspectives of further development of the considered approach are discussed 
in Sect. Conclusion. 
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2 General Definitions 



2.1 Variables 

We consider the reaction where initial particles a and b create 3 particles in the final state. 
We use the notations for momenta of particles and for the set of basic variables which are 
very close to that of the book 0. They are shown in the diagram 



Qi 




Q3 



The choice of labels a and b is made for convenience of classification scheme of the 
forthcoming integration. Momentum transfer t a (t b ) inherits label of the incoming particle. 
For the given particle a ( b ) of the diagram (1) there is the pair of associated variables 
s, t a (s, tb) (the energy variable s is assumed to be fixed in the course of integration) and 
there exists only one variable among the two-particle energies of the final state which is 
nonadjacent to the pair s, t a (s,t b ) — it gets the same label a ( b ). 

Here, we only list the definitions of invariant variables s, s a , s b , t a , t b in terms of 4- 
momenta p a , pb, q%, q 2 , <?3 and masses m a , m b , m 1 , m 2 , m 3 of particles: 



s = 


(Pa+Pb) 2 


= m l 


+ m l 


+ 2p a 


■Pb , 


Sa = 


(<?! + <?2) 2 


= ml 


+ mj 


+ 2 qi ■ 


q2 , 


Sb = 


(92 + 9s) 2 


= ml 


+ mj 


+ 2q 2 - 


qs , 


ta = 


(Pa - qs) 2 


= m l 


+ m t 


-2p a 




tb = 


(Pb - qi) 2 


= m 2 b 


+ m\ 


~2p b - 


qi ■ 



The most general case of particle masses will be considered: all masses are different 
and nonzero — otherwise considerable simplifications are known to take place. 

Every 4-momentum of a particle can enter Lorentz-invariant expression only via scalar 
products with another momenta (we consider the amplitude of the unpolarized experi- 
ment, so there are no other 4- vectors for constructing invariants). 

Five momenta can form 10 scalar products provided the mass-shell conditions 

2 22 22 22 22 2 /o\ 

Pa = m a) p b =m b , q 1 = m 1 , q 2 = m 2 , q 3 = m 3 (3) 

are fulfilled for all external particles of the considered reaction. It is to be noted that in a 
given experiment a specific variable, for example, like (p b — p a ) ■ (qi — q 3 ), might undergo 
investigation. In principle, there is combinatorically large amount of invariants formed 
by linear combinations of scalar products of momenta of 5 particles. Here, we restrict 
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ourselves to the set of invariants shown in the diagram (1). We note only three features 
of the quoted diagram. 

First, these very variables, being invariant masses of pairs of external particles, enter 
the pole contributions of 2-particle resonances. Therefore, these variables are common to 
almost every theoretical analysis. 

Second, the most important property of the above variables is their independence. 
This means that any Lorentz- invariant function of five 4- momenta of the diagram (1) can 
be written in terms of these 5 invariant variables only. Avoiding lengthy proof of the 
independence let us simply give the explicit expressions of all 10 products in terms of 
variables s, s a ,t a , s b ,t b : 

Pb-Pa = (s - m 2 b - m 2 b )/2 , 
Pb-qi = -{h - ml -ml)/2 , 

91- 92 = (s a -ml- m 2 2 )/2 , (4) 

92- 93 = (s b -ml -m\)/2 , 
Pa -93 = -(ta-ml-ml)/2; 

/ 2 , 2 2,2 2 

Pb ■ 92 = [m 1 + m 2 - m 3 + m b - m b 

+2q 1 ■ q 2 - 2p b ■ q 1 + 2p a ■ q 3 )/2 , 
Pb ■ 93 = (—ml — ml + ml + ml + ml 

-2q x ■ q 2 - 2p a ■ q 3 + 2p b ■ p a )/2 , 

/ 2 2 2,2,2 

p a - qi = [m 1 - m 2 - m 3 + m b + m b 

-2q 2 ■ q 3 - 2p b ■ q 1 + 2p b ■ p a )/2 , (5) 

Pa ■ 92 = (- m l + mj + m 3 — m l + m l 

+2q 2 ■ q 3 + 2p b ■ q 1 - 2p a ■ q 3 )/2 , 

/ 2 2 2,2,2 

93 ' 9i = \~ m i ~ m 2 — m 3 + m b + m b 

-2q 2 ■ q 3 - 2q x ■ q 2 + 2p b ■ p a )/2 . 

Here, 5 relations of the first group are simple inversion of definitions @; the relations of 
the second group express the rest scalar products in terms of ones of the first group. 

Third, the planar character of the diagram (1) makes it convenient to introduce the 
notion of adjacent and nonadjacent pairs of variables. In the forthcoming analysis one will 
find a principal difference of the geometry of the 2-dimensional projections of the phase 
space for pairs of adjacent and nonadjacent variables. In the amplitude analysis the 
difference is displayed by the fact that variables of nonadjacent pair can simultaneously 
enter the two-pole contribution as invariant masses of two-particle resonances whereas 
variables of adjacent pair never can meet together in double-pole term. 

To conclude the discussion of the choice of variables one must notice the following: 

1. Along with the set of variables (1) any other set which can be defined in terms of 
a planar diagram will also present the set of independent variables. 

2. Performing substitutions of particles 

(Pb,Pa,Pl,P2,Ps) -» {PiaiPi b iPixiPivPiz) 

with i a , i b , ii,i 2 , i 3 being any transposition of the set {a, b, 1, 2, 3} it is possible to use all 
the results of the analysis of a particular case in any other planar setting. This is provided 
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by invariance properties of kinematical functions determining the geometry of the phase 
space. These properties are discussed in details in the book @. 

3. However, one can not exploit the full 5! = 120 transpositions when dealing with 
specific process, for example, 

7P — > n + n°n . (6) 

Indeed, permutations of particles belonging to initial and final states, for example, 7 <-> 
7T° or p <-> n, introduces variables for another physical processes rather than provides 
description of the considered process in terms of another set of variables. 
Therefore, the principal number of different sets is 2! x 3! = 12. 

4. There are many other sets of invariant variables which, being strongly nonplanar, 
nevertheless admit the same treatment (an example of such a set is: {s, s a ,t a ,t x = (p& — 
Pa) ■ (<?i — 0.2)1 t y = (pb + p a ) 4 (qi — Q2)}) ■ Here, our discussion of the principal problem of 
integrating over phase space will be restricted to the case of variables described by planar 
diagram (1). 



2.2 Phase space integral in invariant variables 

The end point of theoretical analysis of the reaction (1) is the cross section a (a + b — > 
1 + 2 + 3) which might be written in the form 

' = "•//II pfe 2 *^ +p °- qi ~ 92 " 93)|M|2 (7) 

in the case of normalization convention adopted, for example, in the book |J. Here, 
a c = (he) 2 = 0.38937966(23) [GeV 2 mbarn] is the conversion constant, / — statistical 
factor (equal to product of l/n a ! over subsets of identical particles) and 



4J = AJ(p b -p a ) 2 -m 2 a m\ = 2JX(s,ml,m 



-,2\ 



\((p + q)\p\q 2 ) = -A 



p 2 p ■ q 
p ■ q q 2 



(9) 



stands for normalization of initial state. 

In the papers on particle physics different factors in the denominator of the integrand 
(0) like (27r) 3 or even 



2g i0 = 2Vmf + |q/ (10) 

sometimes are used to be hidden into the matrix element by normalization convention for 
particle states. Therefore, the definition of the empty phase space R3 = i?3(|M| 2 = 1) is 
usually based on the common part of (|7|): 

/n^V+Pa-E^) (11) 

J j=l ^jO j 

(energy, or, the same, s variable is assumed to be fixed in (0) and (|i~l"D). 

The differential cross sections are always experimentally known in terms of distribu- 
tions. Thus, the integrations in equations (^) and (|TTD are assumed to be performed over 
all region of allowed momenta in the case of total cross section or over subdomain (slice), 
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cut out by bin bounds in the case of distribution. The dimension of a distribution of the 
discussed 2 — >• 3 process might be 1, 2, 3 or 4. 

This directly follows from the counting of integration variables in expression (|TTD: 
there are three particles in the final state; their 3-momentum components (9 = 3x3 
in total ) are integrated while being restricted by 4-momentum-conservation conditions 
expressed by <5 4 - function — 5 = 9 — 4 degrees of freedom remain. 

Here, we do not consider polarization measurements. Therefore, matrix element |M| 2 
has no dependence on the angle in the plane orthogonal to the beam axis in the laboratory 
frame where the target particle is at rest. Then integration over this variable is easily 
performed providing 2ir factor. 

Hence, maximal dimension of a nontrivial distribution equals 4. The reduction of the 
9-dimensional integral (^) ((|Tl"D) to 4-dimensional form in terms of some sets of variables 
of the most interest is discussed in details in the book In particular, for the case of 
invariant variables (|2|) the result reads: 

° = A/" 7 , 2 W9^ 3(|M| 2 ); R 3 (\M\ 2 ) = 7= = = 2 / 3 (|M| 2 ) ; 
2^A(s,m 2 ,m 2 ) i^f 4 ^A(s, m 2 , m 2 ) 

r 3 (|M| 2 ) = J ds a dtJs b dt b ^0^-\M\ 2 . (12) 

Here, A4 is the Gram determinant of any four independent momenta, say p b ,p a , Qi-,0.2] its 
explicit form in terms of scalar variables (|2|) is given in the Appendix. In what follows we 
omit subscript 3 in (|12]) , (|7]) pointing to the number of particles in the final state. 

In practice, the overwhelming majority of experimental information on rare pro- 
cesses at intermediate energies is represented in the form of total cross sections and 
1-dimensional distributions. The reason might be illustrated by n~p — > 7r~7r + n expe- 



riments []T0| : 1023 full-kinematics events constitute solid ground for total cross section; 
10-14 bins of 1-dimensional distribution have good filling with averaged number of 60-100 
events per bin; the filling of 8 x 8 bins of 2-dimensional distribution is satisfactory (15 
events per bin in average) while filling of 6 x 6 x 6 of 3-dimensional ones and 4x4x4x4 
4-dimensional bins is poor. 

The increase of statistics of the contemporary experiments [1]J by a factor of 10 consid- 



erably improves accuracy of results of total cross sections and 1-dimensional distributions. 
When dealing with 4-dimensional distributions one anyway has to make difficult choice 
between poor binning or insufficient filling of bins. 

Therefore, the minimal problem which solution provides maximal effect is the problem 
of integration over bins of 1-dimensional distributions (an overall integration for total 
cross sections is then solved by simple summation). 

There is well known property (see, for example, @) of the nonadjacent pairs (namely, 
(s a ,t a ) and (s b ,t b )) of variables: projections of the phase space onto a plane of any such 
pair admits relatively simple description. This makes possible to modify considerations in 
such a way that integration over 1-dimensional bins might be in fact realized by adding 
up results of 2-dimensional distribution, known as Chew-Low plot: 



/ 



r A (\M\ 2 ;j,k) = 

ds a dt a ds b dh^^\M\ 2 X (s a] si~\ 4) X (t a ; t k ~\ t k a ) ; (13) 
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x(x; xi,x 2 ) = Q(x - xi)Q(x 2 - x) , 



where s{, t^ 1 , t* are bounds of bin (j, k) in the quoted variables. 

Owing to the symmetry of the phase space in the pairs of variables (s a ,t a ) <-> (s b ,t b ) 
which is induced by transpositions of particles: (a «-> b), (1 <-> 3), it is sufficient to perform 
analysis only for the case (|13"D. To treat 1-dimensional bins in variable of s- (or t— ) type 
one needs to rearrange particles in such a way that the variable in question becomes s a 
(or t a ) variable of the diagram (1) and use trivial binning in the accompanying variable 
of nonadjacent pair — t a (or s a ). 
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3 Geometry of Phase Space 



The present section deals with the equations determining the boundary of the phase space 
over which the integration in eqs. (|?J), flT2p, flT3| ) has to be performed. The main goal 
is to display the origin of expressions for limits of the quoted above integrals when the 
latter are written in an appropriate successive (iterated) form. This requires treating 
unwieldy formulae when they are explicitly expanded. Therefore, we collect the most 
cumbersome final expressions in Appendix, devoting the discussion of the present section 
to the principal steps of analysis of the phase space geometry. 

3.1 General properties of phase space of two nonadjacent vari- 
ables 

The conditions fixing the 3-particle phase space (physical region) of the process (1) might 
be written in terms of inequalities for Gram determinants A n (p 1? ...,p n ) built of all possible 
sets of n momenta of particles: 

A 1 (p 1 )>0, A 2 (p 1 ,p 2 )<0; (14) 
A 3 (pi,p 2 ,p 3 ) >0 ; (15) 
A 4 (pi,P2,P3,P4) < ; (16) 

A 5 (pi,P2,P3,P4,P5) = . (17) 

Here, {p^ stands for any subset of momenta {p a ,Pb, Qi, Q2, Q3} entering diagram (1). 

The thorough analysis of how these conditions arise and should be treated as well as 
references to original papers might be found in J7J. We only remind that 

• left hand sides of conditions (|14D , (|15|) , (|16"D , (|T7|) are invariant functions of momenta 
of particles; substitution of expressions (§), (||) for scalar products of momenta shows that 
eqs. (|14"D, (|16"D, (|17D are conditions for s a ,t a ,Sb,tb at fixed s; 

• conditions fll~4] ) are in fact statements that all pj describe physical particles with 
nonnegative mass squared (Ai(py) = > 0) and positive energy (p j > 0); 

• in the given pattern of iterated integration only few of conditions (pl|) are necessary 
for due treatment (see, for example, ||); 

• the choice of momenta in the determinant A4(pi,p2,P3,P4) has no influence on its 
value — it is the unique function of independent invariant variables; 

• eq. (|17]) expresses the simple fact that in the 4-dimensional Minkowsky space there 
cannot be more than 4 linearly independent vectors. 

The condition (|T6| ) is the natural starting point of analysis. The left hand side of this 
very condition appears in the form fllTl), (fl3|) of the phase space integral when it is rewrit- 
ten in terms of invariant variables. The explicit expression of Z) 4 = — A 4 (g 2 , q3,Pb,Pa) is 
given in the Appendix in terms of expansions 

U 4 — 2^ a a s ,a t ,f3 3 ,l3 t S a t a S b T b 

a a ,a t ,f3 s ,f3t 

= EW?#- (18) 

It is reasonable to consider _D 4 as a function of such subset of variables in which it takes 
the simplest form. Being generally the form of the fourth degree the above expression is 
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only quadratic in any pair of nonadjacent variables. The total list of such pairs among 
variables of diagram (1) is: (s a ,t a ), (s&, t&) and (t a ,tb). (Here, we do not list pairs formed 
with energy s which is fixed during integration.) 

The pair (i a ,t&) can not provide universal treatment since both momentum transfers 
enter the pair and only the energy variables (s a , s&) left. On the other side the choice of 
any pair (s a , t a ) or (s&, U) can provide covering of all cases of 1-dimensional distributions. 

Therefore, assuming that s a and t a acquire some fixed values from the allowed domain 
- Chew-Low plot (the explicit description will be given a little later) — it is convenient 
to write D4 in the matrix form 




d 4 =[ - j •/,.( ) + /,'.( )+/, 1)n . ,19, 

where 

N • (20) 




Let us list the most important properties of the form ([19]) : 
a) determinant of the matrix b 




Det b = s a D 3a /16 , (21) 



where explicit expression for determinant D 3a = A 3 (q 3 ,p b ,p a ) is given by eq. ( |A.11| ) of 



Appendix, is nonnegative in the physical region since two-particle energy s a is positive 
and for the Gram determinant D 3a condition ([15]) is valid; 

b) trace of the matrix b is nonpositive — analyzing the diagonal elements &02 > &20 

of 

b provided by eqs. ( |A.4|) one can recognize A-function expressions for combinations of 



momenta for which the conditions (|14D are fulfilled; 



c) when D4 is transformed to the centered form ( |A.8| ) the free term b c of the latter 



K = 600 - \b T ■ b~ l ■ b = D 2a D 3a /s a . (22) 

is nonnegative in the physical domain of s a , t a variables — this is direct consequence of 
conditions (|TJ), (0) applied to the RHS of eq. (f^); 

d) whenever determinant \b\ becomes zero b c vanishes also; the ratio 6 C /Det b remains 
finite in the physical domain — by comparing ( P2"D and ( A.12j ) one has 



6 c /(Det b) = 16D 2a /s 2 a . (23) 

These properties imply that condition (|16|) determines an ellipse in the (U, Sb) plane. 
Before proceeding with its analysis let us recall that t a and s a variables are assumed to 
be fixed. The location of the (4, Sb) ellipse and orientation of its principal axes strongly 
depends on the latter variables: this is demonstrated by the Fig. 1 where families of 
ellipses are drawn in the (tb, Sb) plane for different values of t a , s a from the allowed 
domain. Horizontal sequences of elliptic curves are obtained for a fixed value of s a - 
from eq. ( |A.7| ) given in the Appendix one can see that si coordinate of the center of the 
ellipse does not depend on t a variable. 
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3.2 Boundaries of phase space for nonadjacent variables 

The standard way to derive the bounds for the (£&, Sf,) variables is to solve first the equation 
-D4 = for a single- variable form ( |A.13| ) (( |A.14| )) which is scanned in the Appendix: 



= ; (24) 



,,. ( K -b tl ±y/b~ t ) 



c(^) = v v 1 ■ (25) 

Here and in what follows we use the pattern of subscript mnemonics (Left-Right), (Down- 
Up), (Ground-High) to mark the solution. Small letters will be used for boundary func- 
tions (like (I, r) in eqs. (0), (p5])) whereas capital ones — for absolute bounds of variables 
in the considered plot. 

The above solutions by virtue of relations ( |A.15| ), ( [A.16Q (see subsect. 2.2 of Appendix) 
are expressed in terms of coefficients bp Bt p t and discriminants b s ,b t . The latter are also 
written out explicitly in eqs. (|A~T9|) , ( [QUI) and ( |A~2TD , (|A~22D of Appendix. Note that 



the plus sign at square root in ( P5|) ( (p4|) ) corresponds to the left bound of t& (sj,) because 
the denominator 6 i2 = &02 (b S 2 = &20) is already stated to be negative. 

Being the product of factors D 3a and D 3as (D 3a and D 3at ) discriminant b s {b t ) vanish, 
first, when D 3a = — this condition determines the absolute bounds of t a , s a pair of 
variables; it will be discussed a little later — and, second, when D 3as = (D 3at = 0). 
The solutions of the latter provide bounds t^, tjf (s^, sf) of interval spanned by 4 (s^) 
variable at given values of t a and s a : 

.l,r {—bsti ± v^t) , s 

h = — (26) 

l,r (—btsi =1= \fb~Z) , s 

Sb - 2^ • (27) 
The explicit expressions for expansion coefficients b s t/3 and bt s /3 of D 3as and D 3at as well as 
for discriminants bt s and b st of these forms are given in Appendix. The important issue is 
the presence of the factor L^a in both discriminants (cf. eqs. ( |A.31| ), (|A.32| ) and (|A.43|) ). 



Analysis of the above discriminants and the properties a)-c) of the form (19) of the 
previous subsection shows that all the necessary and sufficient conditions for existing a 
nondegenerate domain in (tt>, Sf>) plot are collected in the requirement that free term (b c ) is 
positive. Because of restrictions (|T4]), (p~5f ) the only unambiguous splitting of the product 

(0)* 

D 2a > 0; D 3a > 0. 

These conditions determine the Chew-Low plot in s a ,t a variables. Its boundaries cor- 
respond to curves in (s a ,t a ) plane where LHS's of the above inequalities vanish. Because 
of factorization of the expression ( |A.43| ) for Gram determinant D 2a it is easy to see that 
one of the boundary curve is simply the straight line 

s a -(m l +m 2 ) 2 = (28) 

and the other line provided by D 2a 

s a - (mi - m 2 f = (29) 
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is well out of phase space for the square of two-particle energy. 

Here it is where an asymmetry between the momentum transfer t a and energy variable 
s a appears. The lower (left) absolute bound s° for the latter does not depend on initial 
energy s: 

s° a = {m l + m 2 f . (30) 

Another boundary curve in (s a , t a ) plane is provided by vanishing of the D 3a factor 
of b c . This factor, being the second order form of the s a ,t a variables, is presented in 
Appendix in the manner similar to that of D± case. Now, the quadratic form 




D 3a =[ : a ) r )-i- ,s J •( ; j . ,i nn . (3D 

where 

(32) 



are determined by coefficients provided by eqs. ( |A.34| ) of Appendix, is of hyperbolic 



type. The conclusion is not so difficult to arrive to, basing on the properties of this form 
presented in Appendix. In particular, the roots in a variable (provided the accompanying 
variable is fixed) are 

n, h( . v (-(A 11 t a + A 10 )±y/7Q 



2A 20 
A 

2i~ 02 



eM = { -(A u s a + A m) ±VT t) _ 



where discriminants A,, A f are 



A s = (A u t a + A 10 ) 2 - AA 20 (A 02 t 2 a + A 01 t a + A 00 ) (35) 
= D 2c D 2at = D 2c (t a - (m a + m 3 ) 2 )(t a - (m a - m 3 ) 2 )/4 ; 



A t = (A u s a + A 01 ) 2 - AA 02 (A 20 s 2 a + A w s a + A 00 ) (36) 
= D 2c D 2as = D 2c (s a - (Vi + m 3 ) 2 )(s a - (y/s - m 3 ) 2 )/4 . 

This helps to fix critical points t^ ,H , s^ ,H of the considered variables as roots of the 
above discriminants: 

tf H = (m a =F m 3 ) 2 ; (37) 
sf H — {y/s T m 3 ) 2 . (38) 

Unlike the previous case of second-order form in Sb, tb variables, here, the intervals 
between roots are nonphysical since D 3a < there. Regions t a > t% and s a > are 
cross regions of another processes: a + 3— > b + 1 + 2 and a + 6 + 3^1 + 2 respectively. 
It is easy to see that the absolute range of variation of s a variable is fixed in an unique 
way by eqs. ( p0|) and (|38|): 

(*? = s° a ) <s a < (s u a = s f ) . (39) 
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For a given s a bounds tf u (s a ) of t a interval are those given by eq. (|34| ) 

£"(*«) = ti\ Sa ) . (40) 

The absolute bounds of momentum transfer t a are provided in a different way. It is 
simple to find that the lowest possible value of t a is given by the lowest intersection point 
of hyperbola D 3a = with the line s a = 

ta = «(«!) ■ (41) 

For determination of the upper bound it is crucial whether the common point 

si = {s(m a - m 3 ) + m 3 (m 3 m a + m 2 b - m 2 a )}/m a (42) 

of the tangent line t a = and the hyperbola belongs to the physical region (|39|) of s a 
variable or is located below. In the former case the value is never attained. The 
discussed condition depends on the particular relations of the particle masses (and energy 
region); the relevant classification might be found in the paper ||. For the purpose of 
calculations it is sufficient simply to compare with s° a . Then the upper bound t% is 
given by eqs. 

if s T a < s° a then t u a = t h a (s° a ) (43) 
else t u a = t G a . (44) 
The behavior of the s a bounds s^ u at given t a follows the similar scheme: 

s u a (ta) = s h a (t a ) ; (45) 

if t a < t h a (s° a ) then s d a (t a ) = s° a (46) 

else 4(t a ) = si(t a ) . (47) 

The variety of labels is dictated by the need of avoiding overlappings when reproducing 
for the case of (s&, tb)-plot. The symmetry 

a <-> b , 1 <-> 3 (48) 

of diagram (1) allows to use all the formulae of the current subsect. with the substitution 
( fl8D made for all indices and labels containing {a, b, 1, 2, 3}. 

To resume this Sect, one should note that the symmetry in two-particle energy and 
momentum transfer variables which is perfect for s 6 , t b at fixed s a , t a (eqs. (|24]), (p5|), (p7|), 
(p6|) following from the symmetric representation (|i~9"D of boundary function D4) is broken 
in the resulting description of Chew-Low plot: ([39]), (|40|) versus (^3|-4Tf|). This must be 
carefully implemented into algorithms performing integration over phase space. 
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4 Gauss Integration 



4.1 Iterated form of phase space integral 



The analysis of the phase space geometry provided by the previous section makes it 
possible to write the phase space integral ([12]) ( fljjD ) in a number of iterated forms. 
There is the natural splitting of 4-dimensional integral into internal integral over s b , h 
variables 

r b (\M\ 2 ) - 



ds b dt b %^\M\ 2 
and the integral over the rest ones, the final form of (|13|) being 



r A (\M\ 2 ;j,k) 



ds a dt a %J^-x(sa, s{~\ si) X (t a , t k a -\ t k )r b (\M\ 2 ) 



V^4 



(49) 



(50) 



Basing on the calculated boundary functions 
might be written in two iterated forms, namely: 



(H), (|3), (H) the internal integral 



and 



n{\M\ 2 ) 



n{\Mf 



w e 
dt b 



-A, 



I Ml 



f , r>u*») , e(-A 4 ) IWI2 

dt b / ds b ) \M\ 



(51) 



(52) 



i s [{t b ) V-A 4 

There are no reasons to prefer one form or another basing on pure geometric arguments 
- the s b , h domain was stated to be quite symmetric. However, the integrated amplitude 
M might have different behavior in the discussed variables. For example, the amplitude 
might have poles in two-particle energy Sb (shifted to complex plane from the real axe) 
and cuts whereas there should be no such reason of rapid variation with momentum trans- 
fer tb in the integration domain. Since there is maximal step limit among the parameters 
terminating the calculations of integrating procedures it is reasonable to perform the inte- 
gration over "smooth" variable tb first — otherwise all the rest integrations are expected 
to be performed at maximal step number as well. 



As to external integral (50) the difference in the description of (s a ,£ a )-plot in the 
variables which was discussed at the end of previous section makes it more convenient to 
choose the order 



r A (\M\ 2 ;J,k) 



ds n I dt n ^^ X (s a ;si-\si)x(t a ;t k a ~\t k a )r b (\M\ 2 ) 



to simplify the program logic. The order 



(53) 



r A (\M\ 2 ]3 ,k) = 
(Hi! I i " ( ' a) rfSa ^Z^ x(Sa;sr i, si)x( t a; ^-i ; ^) rfc( |M| 2 ) 



(54) 



which directly provide 1-dimensional distribution in momentum transfer t a is less attrac- 
tive from the point of view of amplitude behavior and, what is more important, it causes 
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considerable complications of the program logic which must trace all details of conditions 
(^3j-^7|). Because of these complications the program for calculation of 1-dimensional t a 
distribution based on representation fl54"|) is found to be equivalent to the program for 
processing 2-dimensional distribution according to eq. ([53]). 



When we have written the explicit form fl53|), (51) of the phase space integral (|T3|), the 
problem of calculation of all 1-dimensional distributions in two-particle invariant variables 
is, in principle, solved by an universal procedure which is designed in a straightforward 
way according to the discussed formulae (and boundary eqs. of the previous section 
written in terms of quantities explicitly given in the Appendix). Even with the use of 
integration routines based on Simpson or other standard algorithm the procedure is much 
more effective when compared to the Monte Carlo based analog. 

In the rest part of the current section we shall consider the implementation of the Gauss 
integration method for the integral in question. The possibility of further improvement 
of characteristics of the program comes from the following observations: 

1) the presence of denominator in eq. (^) which vanish at the boundary of the phase 
space makes it evident that the standard integrating routines waste almost all the time 
processing this (integrable) singularity (one can avoid the difficulty by a change of vari- 
ables. This results in minor complications of transfer of arguments from the integrating 
procedure to the user amplitude. It is the loss of universality for further development 
of algorithms processing 3- and 4- dimensional distributions which makes this way of 
correction less attractive); 

2) the integrated matrix element \M\ 2 is usually smooth enough (at least, piecewise) 
function of the variables and it might be well approximated by a polynomial. 

Therefore, calculating the phase space integral in question, one deals with a classical 
case for which the Gauss integration method proved its unprecedented effectiveness. 

Let us note that there might be different approaches of implementation of Gauss 
method to the considered integral. One can try to find the set of orthogonal functions or 
polynomials of all four variables appropriate for the 4-dimensional integral in question. 
This is very interesting and yet unsolved problem. Its solution should help much especially 
in the case when the integration must be performed over all the phase space (for example, 
when analyzing unitarity relation for the considered amplitude). However, when the 
integration is to be performed only over a part (bin) of the phase space or over a sequence 
of bins one easily finds that the set of functions providing calculations by the Gauss 
method must be distinct for every bin. This makes impossible to use once precalculated 
roots and corresponding weights of the polynomial system of the method. 

Here, we shall follow the way which can allow, in principle, to deal with all kinds of 
distributions in invariant variables: 1-, 2-, 3- and 4-dimensional. The separate treatment 
of all four integrals will be found to require only two types of orthogonal polynomials in 
every variable to be considered in the most general case. 
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4.2 Standard singularities of iterated integrals 

To find the polynomial basis which is most suitable for the integral in question one must 
analyse the properties of all four iterated integrals in expressions ([)]]), 

Since it is the phase space induced specifics of the considered integral which is of 
interest we can take the amplitude in the expressions (]5]J) , (|53]) as smooth as we like, say, 
polynomial. 

The singularity structure of the most interior integral over t b in eq. (|5lD is evident. 
The theta-function cuts the entire interval between the roots (t l b , i£) of _D 4 in this variable. 
In the case of a polynomial matrix element \M\ 2 this integral is an elementary one and it 
might be calculated analytically for every monomial term 



V=b^H J(t b -t{)(tl-t b ) 



by recurrent relation 



(55) 



2n - 1 b n n-lb m 

"4-1 T—In-2 ■ (56) 



n 2b t2 n u t2 



This relation uniquely defines I n in terms of coefficients b t o, b t ±, b t2 of the single-variable 
form (|A.14|) of _D 4 . The values of the two starting members of the recurrent sequence 



— ; J i = -7n~ J o ( 57 ) 



7.2 



show that integrand of eq. (0), being polynomial in s b , remains polynomial after first 
integration since the quantity bt 2 , entering both the square root and denominator of 
expressions (f)6|), (|57]) , does not depend on s b - 

This leads to the important conclusions, namely: 

1. The natural weight function for the integral over t b is given by 



Hb{t b 



\l{h-t l b ){ti-t b ) 



-i 



(5f 



2. There is no nontrivial weight function for the next integration over s b . 

It is to be noted that the above conclusions are determined by the integration order 
chosen. The inverted conclusions will be made if one chooses the order of eq. ([52]): it is 
the integral over s b which acquires the weight function, originating from square root of 
D4, the second integral in t b being free from nontrivial weight function. 

In the course of integration via relations (|55]), ( |5"7| ) the irrationality 



H2 



(59) 



appears. It depends on s a , t a variables and might be imagined to be important for analysis 
of the subsequent integrations. In fact, it must be modified by the following integration 
over s b and the final answer appears to be more symmetric in terms of components of the 
elliptic form (|19|). 

Let us shift the s b , t b variables to make the answer less immense. The coordinates s^, 
tl of the ellipse center do depend on the s a ,t a variables. Fortunately, they were found 
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to be smooth in the physical domain (see eqs. ( |A.7| )). The standard monomial in both 
variables is then 

imi 2 = M nm = ( Sb - sxiu - nr . (60) 

By means of elementary calculations the details of which it is reasonable to omit, the 
integral 

r nm = r b {M nm ) (61) 
can be brought to the form (provided n + m is even; otherwise integral r nm is zero): 



2tt 



. n— m . . 

6 02 )— x (62) 



(n + m + 1)!! V Det b \Vet b) 

[?] 

x J2 Cf n {n + m-2l- 1)!!(2Z - l)!!(Det S)' [b n /2} m ~ 2i 
i=i 

or, changing the order of integration, to the equivalent form: 



2tt b c b c 



n-\-m 
2 



(n + m + 1)!! V Det b \BetbJ 

[f] 

x Y,C 2 n l (n + m-2l- 1)!!(2/ - l)!!(Det Vf [6 u /2] 



. m — n 

2o)— x (63) 



At first glance the quantity 602 (620) appears in the denominator of eq. (|62])( (p3D ) when 
m > n (n > m). Then eq. (|B3|) ( (|5^D ) helps to avoid the algebraic proof that the quantity 
in the denominator is exactly cancelled by the factor from the sum. 

Taking into account that coefficients b a /3, entering the answer, are (second order) 
polynomials of the s a , t a variables, one easily continues analysis of properties of the 
iterated integrals: 

3. The true irrationality generated by integration over s&, t& variables is 



b c V [s a ~ (mi + m 2 ) 2 ] [s a - (m 1 - m 2 ) 2 ] 

2 y_i v — — i — . (64) 



Det b 



4. Polynomial matrix element in t a variable remains polynomial after the considered 
internal integration; there is no nontrivial weight function for the next integration in t a 
variable. 

5. The discussed in the previous section asymmetry of the description of the phase 
space in the variables s a , t a appears to be deeper after the internal integration; the generic 
polynomial of s a variable acquires irrationality with the critical point just at the boundary 
of (s a , t a )-plot 

s a = (mi + m 2 ) 2 = s° a (65) 

multiplied by a rational function of s a . 

6. The position of the only pole of this rational function 

s a = (66) 
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is distant from the physical region, an expansion or polynomial approximation being 
allowed (in the potentially dangerous case of vanishing masses mi, m 2 both the pole and 
the irrationality disappear — see eq. fl64|)). 

Turning to the integration over t a variable one should remind that integration in eq. 
(|53|) is assumed to be performed over domain cut out from the (s a , t a )-plot by rectangle 
defined by bounds s J a , S t„ of bin (j, k). Whenever the rectangle is exactly inside 
the Chew-Low plot integration of any t a polynomial will provide only constant (in s a ) 
contribution. If there is an interval in t a for which an arc of hyperbola D 3a = enters 
the integration domain as a part of boundary the t a integration results in an expression 
containing the term \[A~ t from the boundary function fl34|). The expression (|36|) for dis- 



criminant At shows that it is the absolute bound = (y/s — m 3 ) 2 where the analyticity 
might be lost due to the above square root — another critical point = (y/s + m 3 ) 2 is 

bringing to the 



well outside the physical region. (And again, the condition m 3 
boundary simultaneously kills the square root itself.) 

Depending on the given binning in the s a , t a variables, the absolute bounds 
might be unattainable at all, only one might be attained and, at last, for a specific bin 
of a t a distribution in the case of some energy and particle masses both bounds might 
happen to undergo due counting. The bounds must always be taken into account when 
the total cross section is being calculated. 

Combining with the point 5. of the previous conclusion one can find that the last 
integration in ([53]) must be performed either 

a) with the two-point singular expression 



mi + m 2 J 



m 3 j 



(67) 



or 



or 



b) with 



c) with 



yj[s a - (mi +m 2 ) 2 ] 



(Vs- m 3 ) 2 - s a 



(6* 



(69) 



or 



d) without square root singularity. 

Now, it is easy to realize that cases b) and c) admit a simple change of variables 



s a s a X 



(70) 
(71) 

respectively, which helps to get rid of singularity. Of course, the price is the doubling 
of the effective power of integrand (i.e. the power of approximating polynomial for the 
smooth factor of integrand). This looks quite acceptable. In contrast, the case a) needs 
a nonpolynomial (trigonometric) substitution providing no guarantee for estimated rate 
of convergence. 

Finally, we can state that 
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7. Depending on the given binning the last integration has the weight function given 
by eq. f|67|) , namely, 



or might be processed with the trivial weight function; in the latter case a change of 
variable might be required. 

To resume the current subsection one should reread the discussed above points 1., 
7. and note that there is no principal difference between the weight function fibih) 
of eq. (|58| ) and iia(so)- Indeed, by multiplication of numerator and denumerator of 
integrand by ^a{so) one increases the effective power of the latter by one and converts the 
irrationality to the form of fib- Therefore, only two distinct sets of orthogonal polynomials 
are required for implementation of the Gauss method to the integration over the phase 
space in question, namely, the set, corresponding to the (normalized to the unit interval) 
weight function of eq. (0), and the set, corresponding to the unit weight function. 




(72) 
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5 Principal Features of FORTRAN Code 



5.1 Outlines of the Gauss method; specifics of realization 

To proceed with implementation of results of the previous analysis let us briefly recall 
the basic formulae of the Gauss method ||12|| . More details might be found in almost any 



book on numeric integrating. For example, see [13 



Suppose we have to calculate integral of the type 

f(x)fj,(x)dx , (73) 



where f(x) is a smooth enough function and /x — positive weight function which might 
have integrable singularities at the points ±1. Such weight function uniquely defines 
the set of polynomials {<p n }, (n > 0) orthogonal to each other with this weight function 
provided a normalization convention is adopted. Then the Gauss formula reads 

f f(xMx)dx ~ £ 4 n V(4 n) ) = Gn . (74) 

Here, the weights {w^} and points {xjj. £ (—1,1)} are assumed to be chosen in such 
a way that relation fl74|) turns into equality for any / being a linear combination of the 
first n + 1 functions m from the system {0 m }. In this case x^ should be the roots of 
the polynomial (f) n . 

The function / being a smooth function possesses rapidly decreasing Fourier coeffici- 
ents c n : 

f( x ) = c n<Pn(x) ; (75) 

n>0 
»1 

f(x)(f> m fi(x)dx . (76) 

i 

The convergence rate of the integral sums {G n } depends on that of {c n } and, consequen- 
tly, on the choice of the system {0 n } — this opens the possibility to hide singularities of 
an overall integrand into the weight function n(x). 

According to the conclusion of the previous section the two cases of weight functions 
are peculiar to the iterated integrals over phase space, namely, 

H(x) = (1 - x 2 )- 1 ' 2 (77) 

and 

IM>(x) = 1 . (78) 
In the first case the orthogonal functions are the Chebyshev polynomials 

T n (x) = cos (n arccosx) , (79) 

for which there are exact expressions for roots and weights: 

M 2& — 1 
xi = cos Ti ; (80) 
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w 



(«) 



7T 



n 



(81) 



the quadrature formula ( 74|) being 

n / 

f{x){l - x 2 )~ 1/2 dx ~ £ -/ cos 



i 

-l 



fc=i 



2fc 



7T- 



2n 



(82) 



The unit weight function defines the set of the well-known Legendre polynomials 

1 <f 



Pjx) 



2 n n! cfe r 



1 - x 



2\n 



83) 



In this case the roots X^' of the Legendre polynomial P n for n > 10 are to be found only 
by numerical solution of algebraic equation P n (x) = 0. The corresponding weig hts Wi n) 
are then given by the expression 

-2 



w, 



(») 



2[1 - (X 



Explicit expression for the considered integral is 



f(x)dx 



2£[1 

fe=i 



fc , 



dx 



Pn 



X 



(n) 



f(x 



(n)^ 
fc , 



54) 



The formulae (]82f) , ( j85|) are the basis of the central subroutines of the integrating 
program, namely, CHIN and ZNIN. To avoid problems with recurrent calls in the course of 
processing iterated integrals the subroutines were simply cloned: identical copies CHIN1, 
CHIN4 and ZHIN1, ZHIN2, ZHIN3 were used to be called at the corresponding level of 
integration. 

Whereas the programming of a subroutine of the CHIN type (based on Chebyshev 
polynomials) is straightforward ZHIN subroutines require calculation of roots of Legendre 
polynomials and the corresponding weights in much more complicated manner. Therefore, 
to avoid significant losses of time during integration runs the way of using precalculated 
values was chosen. 

The roots of Legendre polynomials have the following property: any root of every 
subsequent polynomial is located between corresponding roots of the previous Legendre 
polynomial. Hence, the Newton method of determination of roots is the most suitable in 
this case. It operates with function and its derivative. Legendre polynomials Pi and the 
derivatives Di might be easily calculated by the well-known recurrent formulae: 

lP t {x) = a;P,_i(z)(2/ - 1) - P,_ 2 (z)(Z - 1) ; (86) 

D t {x) = (21 - l)fU(z) + A- 2 (x) , (87) 

starting the recurrence from Pq = 1, Pi = x and D — 0, D\ — 1. 

As soon as the roots are found the values of corresponding weights are determined by 
eq. (|8~4| ) in terms of roots and derivatives fl8"7p . 

Since the above calculations are needed to be performed only once the time of run is not 
of much importance but the precision is, because it directly determines the accuracy of the 
integration method. To avoid complicated analysis of the influence of the uncertainties of 
roots and weights on the final answer the latter were calculated at the maximal precision 
allowed by VAX/VMS PASCAL. The present realization contains common block RTWT 
with roots and weights for Legendre polynomials up to 50 order (up to 100 — available); 
in practice, the maximal order 16 appears to be called. 
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5.2 Implementation 

The principle structure of the program is simple: User defined MAIN manipulates with 
User data base to provide necessary parameters, calls integration subroutine WTFF and 
performs appropriate output of results. There must be another User routine which is 
called by WTFF, namely, matrix element USRF. 

The numeric integration over 4-dimensional phase space is performed by subroutine 
WTFF. It returns the value of integral and its absolute difference with the value of previous 
iteration which are the variables of the argument list. The list of physical variables to 
be passed to this procedure includes energy s, particle masses m a , nib, mi, m 2 , m% and 
bounds of 2-dimensional bin s{, t^ -1 , t*. Only the bounds were chosen to be included 
into the list of arguments — other variables are accessible through the COMMON blocks. 
This is because of the following reasons: a) along with the value of s some other totally 
equivalent precalculated quantities Tbeam, -Ebeam, -Pbeam, etc. are necessary also for 
calculation of matrix element provided by User; b) according to the described in the end 
of section 2 universal way of treating bins in one or another variable the above list of 
masses realizes some transposition of particle masses of the User reaction; a manipulation 
is needed to pass true arguments to the amplitude which should not know anything about 
this. 

Apart the above physical variables the subroutine WTFF operates with a number of 
control parameters including terminating value EPS of attained accuracy and 4 limits of 
iterations for all 4 integrals. 

The WTFF algorithm exploits the crucial property of phase space which is expressed by 
the fact that at any given values of variables s a and t a the allowed domain of variables s&, 
tf, is the ellipse. Parameters of the latter are determined by masses of external particles 
and the values of s a and t a — the corresponding functions and subroutines are built 
according to the formulae given in section 3 and in Appendix. 

The integration domain in the plain (s a ,t a ) is the intersection of the given rectangle 

st 1 <s a < si ; (88) 

<ta<t k a 

with the region 

< *a < s u a ; (89) 
t a (s a ) <t a < t h a (s a ) 

(see eqs. ©, ©). 

The properties of the resulting domain determine the integration method (based on 
Legendre or Chebyshev polynomials) to be applied to the s a variable. In the most general 
case we have to split exterior integration into up to 4 processes (intervals) (SINTA, 
SINTD) so that boundaries of the internal integration over t a have smooth dependence 
on s a on each interval. Thus every process performs integration of the smooth function 
FWA, FWB, FWC or FWD (which are specified below) and admits an exact method. 

As it is discussed in subsection 4.2 almost each of the four intervals might contain 
singularities at their ends classified by the cases a), b), c) and d). So, depending on the 
relative values of integration interval bounds and absolute bounds (^), the calculation of 
each of integrals SINTA, SINTD is processed by different algorithms, using procedures 
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CHIN1 or ZHIN1. The algorithms call the corresponding functions FWA, FWB, FWC or 

FWD depending on the case of singularities classified by a), b), c), d). These functions 
realize appropriate changes of variable s a and make this problem hidden for the next 
integration over variable t a . 

The analysis of the previous section (see points 1., 7.) shows that the next two 
integrations over any bin in variable t a and over the complete range (sf , ) of Sb variable 
are performed by Gauss method with Legendre polynomials. The algorithms in question 
are realized by identical subroutines ZHIN2 and ZHIN3. 

According to the quoted above analysis the last integration is processed by subroutine 
CHIN4 which realizes a simple algorithm of Gauss integration with Chebyshev polynomials 
over unnormalized interval (t l b , t r b ) given by eq. fl2~5l) . The call of the User matrix element 
USRF of (s a ,t a , Sb,tt)) arguments arises at this stage. Prior to the call the initialization of 
all ten scalar products (|4]), (Q) is performed. The scalar products (like the particle masses) 
are available via COMMON block. This must simplify the calculation of the considered 
amplitude which is usually derived according to Feynman rules. However, User must 
provide appropriate rearrangement of the scalar products if variables of the diagram (1) 
are a permutation of the User ones. 

In total, the discussed program contains about four dozens of functions and subrou- 
tines handling various kinematic calculations. The detailed description will follow the 
specifics of the phase space geometry discussed in Sect. 3 and collection of formulae given 
in Appendix. 
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6 Conclusion 



In the result of analysis of the geometry of relativistic 3-particle phase space and im- 
plementation of integration algorithm the principal solvability of the integration problem 
appears to be demonstrated. Besides, it becomes evident that the considered case of 2 
— > 3 reaction is likely to be the last one permitting a treatment outside the framework of 
Monte Carlo approach. This is due to combinatorical growth of the number of kinemati- 
cal variables (as free as well as dependent ones) characterizing the reaction (as compared 
with 2^2 case) and the complication of geometry which is displayed by growing alge- 
braic power of definition equations. Given more than 3 particles in the final state other 
methods of analysis, like the artificial intelligence ones, become necessary. 

The integration problem has been solved here for 1 -dimensional distributions in any 
of 10 — 1 = 9 two-particle invariant variables (i.e. without counting initial energy s) of 
the considered reaction on the base of the fast integration algorithm implementing the 
Gauss method. The algorithm needs to process the calculation of matrix element in few 
thousands of points in phase space (compared with tens of millions ones and even more 
in the Monte Carlo routines) and provides the accuracy which will be never attainable by 
Monte Carlo calculations. 

It goes without saying that the solution was possible due to the extensive coverage of 
the 2^3 case by preceding investigations (the reader had been already advised to look 
the books |7], [J for an overview). The same reason might be used for motivation of further 
researches in this field. At least two directions of the researches might be proposed: 1) 
a generalization of the approach for the case of noninvariant variables (like angles and 
3-momenta); 2) its extension to treat 2-, 3- and 4-dimensional distributions. 

While the first direction seems to be closer to the needs of specific experiments the 
second one appears to be more attractive. Indeed, now it is not so difficult to form almost 
any distribution (hence, in invariant variables as well) from full-kinematics data with the 
help of modern tools like GEANT-3. The capability of comparison of different measure- 
ments will be opened then and the phenomenological parameters describing specifics of 
distributions will get an universal status. 

The field of application of the considered integration algorithm might include a large 
variety of processes at intermediate energies like nD — > npn, tcN — > ttttN, jp — > 7nrN, 
etc. One can also try the approach at high energies as well (at least for control of accuracy 
of the ultra-relativistic approximation for phase space). 

7 Acknowledgments 

We thank CHAOS team and DEC for providing resources of ALPHA processors at TRI- 
UMF (Vancouver), INFN (Trieste) and AXP/OSF (Pasadena). AAB is grateful to col- 
leagues of CHAOS team for interest, hospitality and help during his stay at Canada and 
Italy. 

This research was supported in part by RFBR grant N 95-02-05574a. 



24 



References 

[1] R. Brun, et al, "GEANT-3" , CERN preprint DD/EE/84-1 (1987). 

[2] R.A.Arndt, J.M.Ford and L.D. Roper. Phys. Rev., D32 (1985) 1085. 

[3] A. A. Bolokhov, V.V. Vereshagin and S.G.Sherman. "Analytic integration over ttN — > 
ititN reaction phase space". LNPI preprint N 1542, Leningrad, 1989. 

[4] D. Bardin, A. Leike and T. Riemann, "Semi-analytical approach to Higgs production 
at LEP-2". DESY preprint 94-097, 1994. 

[5] G.F. Chew and F.E. Low, Phys. Rev. 113, 1640(1959) 

[6] G.I. Kopylov, "Foundations of resonance kinematics", Nauka, Moscow (1970) (in 
russian). 

[7] E. Byckling and K. Kajantie , "Particle Kinematics", John Wiley and sons, N.Y. 
(1973). 

[8] K. Kajantie and PLindblom, Phys. Rev. 175 (1968) 2203. 

[9] J.F. Donoghue, E. Golowich and B. Holstein: Dynamics of the Standard Model. 
Cambridge University Press, NY, 1992 

[10] T.D.Blokhintseva et al., JETP, 44, 498, 1963. 

[11] G. Smith and M. Kermani, "Experiment 624: The (tt,2tc) reaction, a tool to deter- 
mine scattering lengths and coupling constants" . TRIUMF Annual Report Scientific 
Activities 1994. 

[12] C.F. Gauss, Werke, Bd3, Gott., 1866, S. 163-196. 

[13] A.H. Stroud and D. Secrest, Gaussian Quadrature Formulas, N.Y., 1966. 

[14] R.A. Morrow, J. Math. Phys. 7 844(1966). 

[15] J. Tarski, J. Math. Phys. 1 149(1960). 



25 



8 Appendix: Gram Determinants 



8.1 Two— variable matrix form of D4 

The fourth order Gram determinant -D4 does not depend on the choice of particle mo- 
menta. Let us define it as follows: 

£> 4 = -A 4 (g 2 ,g3,Pb,Pa) 



92 


92 


92 


93 


92 


Pb 


92 


■Pa 


93 


92 


93 


93 


93 


Pb 


93 


■Pa 


Pb 


92 


Pb 


93 


Pb 


Pb 


Pb 


■Pa 


Pa 


92 


Pa 


93 


Pa 


■Pb 


Pa 


■Pa 



(A.l) 



Here, all scalar products depending on Sf, and U variables are underlined (see eqs. (H), 
(H)). This makes evident that is only quadratic in Sb and tb- By virtue of the symmetry 
(a «-> b, 1 «-> 3) the same conclusion is valid for (s a , t a ) pair as well. 
Expression for D± in terms of expansion in s a , t a , Sb, h: 

D 4 = J2 d as , aM s a a s CsK (A.2) 

ot s ,at,/3 s ,/3t 

= (A.3) 

is determined by nonzero coefficients 

doooo = (-(2(((m 2 b + w 2 - s)m 2 3 + m\m\ - + m\s)m\ - m\m\ 
+m\m 2 m 2 a )m\ + 2(((m 2 + s) — m\)m\ — 2m 2 a s)m\m b 1 
+ (m 3 + m a ) 2 (m 3 - m a ) 2 m^ 

+ (m 2 + 2m b m a + m 2 - s)(ml - 2m b m a + m 2 a - 5)777,3 
+mimt))/16 ; 

^0001 = (((2777^ - 777 2 - S)m 2 + 777^ + 777 2 S)777 2 

// 2 2 , \ 2 , 2, 2 2\ 2 

— ((777 fe — 777 a + Sj777 3 + 777 fc S + 777 a S — S J777 2 

+777 3 777 6 — 777 3 777 6 s)/8 J 

^0002 = - s) 2 )/16 ; 

^0010 = (-{{{ml + S)- 27T7 2 - 777 2 )7772 

+ (m 3 + r77 a )(r77 3 - 777 a )777 2 - 777 2 777 2 )t77 2 ) /8 J 

doon = - s ) m &)/8 ; 

^0020 = (— r77.f)/16 ; 

rfoioo = (((ml - 777 2 - s)m 2 - (m 2 b - 2m 2 + s)m 2 3 - m 2 a s)m\ 

,/2, 2\ 4 / 2 1 2 \ 2 , 22\/o 

+ (777 3 + m a )m 1 — [m b + 777 a — sj777 2 s + ■m 3 m b s)/S ; 

^0101 = ( — ((777 3 + s)777 2 — 27772^ + r77 3 S — S 2 ))/8 ; 
^0110 = (-((777 2 + 777^ + m 2 -2s)777 2 

+ (777 2 - 777 2 + s)777 2 + 777 3 777 2 + 777 2 s))/8 ; 

4n = (mg - s)/8 ; 
doi2o = m 2 b /8 ; 
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^0200 = {-{ml - s) 2 )/16 ; 

^0210 = (ml + s)/8; 

^0220 = (— 1)/16 ; 

diooo = {{{^ml-m 2 3 + m 2 a )ml-{m 2 b -m 2 a + s)ml + m 2 3 m 2 b )m 2 a )/8; 

diooi = (-(("is + m l) m i + ( m t + m l ~ 2s)ml 

-{m 2 b -m 2 a - s)m\ + m^s))/8 ; 
c?ioo2 = (mg + s)/8; 
dioio = (( m 3 + m a ){m 3 - m a )m 2 

/ 2 , 2 \ 2 22 2 2\/o 

-{m b + m a - s)m 2 - m 3 m b - m b m a )/8 ; 

dion = {-{ml + ml - 2m 2 a + s))/8 ; 

^1020 = ; 

dnoo = (-(m? - s)m^)/8 ; 

rfnoi = {mj - s)/8 ; 

dmo = (-K-2m^ + m^ + s))/8; 

dun = (~l)/8; 

^1120 = 1/8 ; 

tfeooo = (-m^)/16 ; 

^2001 = m 2 a /8 ; 

G?2002 = ( — 1)/16 ; 

tfeoio = m 2 j8 ; 

<^2on = 1/8 ; 

G?2020 = ( — 1)/16 . 

Coefficients &g S)/ 3 t are built of combinations of the latter ones with corresponding powers 
of {s a ,t a ): 

b 00 = {-{2{{{m 2 b +m 2 a -s)m 2 3 + {m 2 a -t a )m 2 b 
+{s - 2s a + t a )m 2 a -m^ + st a )m\ 
+ {{m 2 a + t a )m 2 b + {s a - 2t a )m 2 a + st a )m 2 3 
+{s + s a )m 2 a t a - m\m 2 b - m\s a - st 2 a )ml 
+2{{{m 2 a + s) - m 2 b )m\m 2 b - ((2s - s a )m 2 a - st a )m 2 b 
+ {s a + t a )m 2 a s - mls a - sH a )m\ 

+ {m\ + 2m 3 m a + m 2 a - t a ){m 2 3 - 2m 3 m a + m 2 a - t a )m\ 
+ {ml + 2m b m a + m 2 a - s)(m\ - 2m b m a + m 2 a - s)m 2 
-2{m 2 a s a + st a )m\m 2 b + m\m\ + m\s 2 a - 2m 2 a ss a t a + s 2 t 2 a ))/lQ ; 
6 i = {{{2m 2 b - m 2 - s - s a - t a )ml + (s - s a )m 2 a + m\ - st a + s a t a )ml 
-{{m 2 b -m 2 a + s)m\ + {s + s a )m 2 a + {s - s a )m 2 b 
+ {s a - 2t a )s - s 2 )m 2 2 

-{{s + s a )m 2 b - {2s a - t a )s + m 2 a s a )m 2 3 - {s - s a )m 2 a s a 
+m\ml + s 2 t a - ss a t a )/8 ; 
b io = {-{{{m 2 b - s a + t a )ml - {m 2 a - t a )m 2 b 
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+ {s a + t a )m 2 a - 2st a + s a t a - t 2 a )m\ 

+ ({m 2 a + s + s a + t a )ml + (m 2 a - s)(s a - t a ) - 2m 2 3 m 2 b - m A b )m 2 2 

+((s a + t a ) - m 2 b )mlm 2 b + {m 2 a s a + st a - 2s a t a )m 2 b 

-(m 2 a s a - st a )(s a - t )))/8 ; 
hi = (-(((s a -t a ) +ml)ml- (s - s a )m 2 b + (s a + t a )s 

~(s a ~ t a )s a - 2m 2 a s a ))/8 ; 
bo2 = (2(s + s a )ml - m\ - s 2 + 2ss a - si)/ 16 ; 

b 20 = (2(s» + t a )ml - ml - s 2 a + 2s a t a - t 2 a )/16 . (A.4) 
When considered as form of (s&,£&) variables -D4 might be written as 

M::) T ' s '0:) + » T 'U') + ^ ( ^ 



where 

b= (b W ) . h= ( b 20 b u /2\ 
Eliminating the linear terms in the form ( |A.5[ ) by shift of variables 



(A.6) 



where 



--r 1 • b (A.7) 



V n J 2 

( {~{s- s a - m\)m\ + (s + s a - m 2 3 )m\ + (s - s a + m§)s )/ (2s ) \ 

V + *a - ™fc) m l + (S a - t a + ^b) TO 2 ~ ( S a ~ *a ~ mg)s a )/(2s a ) J ' 

it is possible to rewrite D4 as the centered form 

D4= u) T - 6 -(?n +6 " <A - 8) 



The free term b r of this form is 



b c = b 00 - ~6 T -b- l -b = D 2a D 3a /s a . (A.9) 

Here, D 2a and D 3a are the Gram determinants 

D 2a = -A 2 (q u q 2 ) = (s a - (m 1 +m 2 ) 2 )(s a - (m 1 -m 2 ) 2 )/4: , (A.10) 

D 3a = A 3 (q 3 ,p b ,p a ) 

= {-s a m 4 a - t a s 2 + t a s(s a - t a ) 

+ {{m 2 a + s + s a + t a )m 2 b + {m 2 a - s){s a - t a ) - m^m 2 (A. 11) 

+((s a + t a )s - s 2 a + s a t a )m 2 a - (m 2 a - t a )(s - s a )m 2 b - m\m 2 b )/^ . 



In its turn the determinant of the matrix b of the quadratic form QA.5 ) which remains the 
same for the centered case (|A.8|) also contains the D 3a factor: 

Det b = s a D 3a /l6 . (A. 12) 
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8.2 Single— variable form of 

When considered as function of a single variable (because of symmetry a <-> b the variable 
is reasonable to choose s b or t b ) _D 4 is determined by expansions 

£>4 = b s2 s 2 b + b sl s b + b s0 , (A.13) 

£>4 = Mft+Mt + fctO , (A. 14) 



where 



and 



b S 2 — b 2 o ; 

6.1 = 6ut 6 + 6io; (A. 15) 

b s o = b 02 tl + b it b + b 00 



bt2 — bo2 ', 

6*i = 6ns 6 + b i ; (A.16) 



ho = b 20 s 2 b + b 10 s b + b 



are defined by expressions ( |A.4| ) for coefficients bp s> p t . 

To calculate the roots of the second order polynomials ([A. 13 ), ( A.14|) one needs in 



fact only coefficients b s2 and b s i (b t2 and b t \) provided the discriminants of ( |A.13j ), (|A.14j) 

b s = b 2 sl -Ab s2 b s0 , (A. 17) 

b t = b 2 n -4b t2 b t0 (A. 18) 

are known. The relevant discussion of Jacobi reduction theorem for symmetric determi- 



nants of arbitrary dimension might be found, for example in |TJ], [HJ. In our particular 
case of the fourth order Gram determinant the results might be verified by direct 
calculation to be the products of Gram determinants 

b s = D 3a D 3as , (A. 19) 

h = D 3a D 3at , (A.20) 
where explicit expression for D 3a is given in (|A.11[) and the ones for D 3as , D 3at are 

D 3as = {{{m 2 b - s a + t a )m 2 2 - (t b -t a )m 2 b + (s a + t a )t b + s a t a -t 2 a )ml 

+ s a + t a )m 2 + (s a - t a )t b - mi)ml (A.21) 
+ (t b - t a )m 2 b s a - (s a - t a )t b s a - m\t a - m 2 m 2 h - t 2 b s a }/A ; 

D 3a t = {(((^3 + s - s a )m 2 2 + + s a + s b )ml- m| - ss b + s a s b )m\ 

+((s - s b )ml + (s a + s b )s - s 2 + s a s b )m 2 2 (A.22) 
-(s - s b )mls a - m\m\ - m\s + ss a s b - s 2 a s b - s a s 2 b )}/A . 

For processing the boundaries of (s b , t b )-plot at given fixed values of s a and t a variables 
the similar treatment of D 3as and D 3at determinants is necessary. Expanding the latter 
in the form 

Dzas = b st2 t\ + bstih + b st Q ; (A.23) 
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where 



D 



3at 



b tS 2S 2 b + b tsl s b + b 



tsO , 



(A.24) 



b s to = {((ml - s a + t a )ml + m 2 b t a + s a t a - t 2 a )m\ 

+ ((s a + t a ) - ml)m\m\ - m\t a - m\m 2 b - m 2 b s a t a )}/A 

b s ti = {{ml + s a - t a )ml - (ml - s a - t a )m\ 
+m 2 b s a - s 2 a + s a t a }/A ; 

b st 2 = -s a /4 ; 



(A.25) 
(A.26) 

(A.27) 



btso 
hsi 
b tS 2 



{((ml + s - s a )m\ + (s + s a )m\ - m\)m\ 

+ (777.3 — s + s a)mls — m\m\ — m\s — m\sSa] /4 ; 

{(ml - s + s a )m\ - (ml- s- s a )m 2 2 

+mls a + ss a - s 2 a }/A ; 



(A.28) 
(A.29) 

(A.30) 



along with coefficients b st2 , b st i, b ts2 , b ts i one needs to know the discriminants of ( |A.23| ) 

(Em- 



= b 



stl 



46,«6 



st2<->st0 



-(2(s a + t a )m 2 b -mt-s 2 a + 2s a t a - t 2 a ) x 



(m 1 + 2m 1 m 2 



m 



s a )(m\ - 2m x m 2 + m\ - s a )/lQ 



(A.31) 



b tsl - 4:b tS 2btso 



-(2(s + s a )ml -m\- s 2 + 2ss a - s 2 a ) x 

(m\ + 2777x7772 + m 2 — s a)( m l — 2mim 2 + 777-2 — s a)/16 



(A.32) 



8.3 Two— variable matrix form of D 



3 a 



It is again convenient to determine the Gram determinant D 3a defined by eq. ( |A.11| ) as 
expansion in its variables s a ,t a : 



D 3a = £ A aa 
Explicit expressions for coefficients are 



Q Ot t -lOlt 

,at s a °a 



(A.33) 



^00 = 


((((m 2 a + s) 


A01 = 


1 1 2 2 
((m b - m a - 


A w = 


((m 2 b + m 2 a - 


A 02 = 




A n = 


(—(777^ — 777 


A 20 = 


(-m 2 a )/4 . 



i 2 b s + 

,2^2 



s 2 )/4 

1 ™2„\ 



s))/A 



(A.34) 
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The matrix A of the quadratic form ( |A.33| ) 



D -=(t:) T - A {Z) +AT { s t :) +A °°' (A - 35) 

where 

has the determinant which up to a constant coincides with the Gram determinant 

D 2c = -A 2 (p a ,p b ) = (s - (m a + m b ) 2 )(s - (m a - m b ) 2 )/4: , (A.37) 

namely, 

A 2 = Det A = -D 2c /1Q . (A.38) 
The free term A c of the corresponding to ( |A.35| ) centered form 

A c = A 00 - ^A T -A- 1 -A (A.39) 

also is expressed via D 2c : 

A c = -m\D 2c = -ml(s - (m a + m b ) 2 )(s - (m« - m b f)/A . (A.40) 

Let us give for completeness the coordinates of geometrical center of the conic section 
determined by the form D^ a : 

(S)--^-(^)- (A ' 41) 

8.4 Second— order Gram determinants 

To control absolute bounds of s b , t b variables one should add to the list of determinants 
D 2c , D 2a , D 2asi D 2at 

D 2c = -A 2 ( Pa ,p b ) (A.42) 

= (ml + 2m b m a + m 2 a - s)(ml - 2m b m a + m 2 a - s)/4 ; 

D 2a = -A 2 (q u q 2 ) (A.43) 

= (m\ + 2m 2 mi + m\ — s a ) (m\ — 2m 2 m 1 + m\ — s tt )/4 ; 

D 2as = -A 2 ( qi + q 2 ,q 3 ) (AAA) 

= (ml + 2m Z \fs + s — s tt )(m| - 2m 3 y/s + s - s tt )/4; 

^ 2ai = -A 2 (p a ,g 3 ) (A.45) 

= (ml + 2m 3 m a + m 2 a - t a ) (ml - 2m 3 m a + m\ - t a ) /4 

already used in previous considerations the analogues of D 2a , D 2as , D 2at obtained by 
transposition of particles 1 «-> 3, a <-> 6. 

All these determinants treated as functions of quantities entering their bodies are in 
fact represented by the same standard kinematical function A ( cf. eq. @) 

A 2 (p,q) = Det( p P q2 q ) = -\\((p + q) 2 ,p 2 ,q 2 ) (A.46) 
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with 



X(x, y, z) = x 2 + y 2 + z 2 - 2xy - 2yz - 2zx (A.47) 
= (x-(v^+v^) 2 )(x-(W-v^) 2 ) 

= (z-(^ + Vy) 2 ){z-(V^-Vy) 2 ) ■ 

The same applies to D 3a , D 3as , D 3at (and, generally speaking, to D 4 as well) — there 
are many forms in which the so called kinematical G-function might be represented and 
many elegant properties of it which then follow. Having in mind quite utilitarian goal we 
display only few of them in the simplest possible terms. For more details reader can look 
into 0. 
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